The Crossover from the Bulk to the Few-Electron limit in Ultrasmall Metallic Grains 
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We study the properties of ultrasmall metallic grains with sizes in the range of 20 up to 400 
electrons. Using a particle-hole version of the DMRG method we compute condensation energies, 
spectroscopic gaps, pairing parameters and particle-hole probabilities of the ground state wave 
function. The results presented in this paper confirm that the bulk superconducting regime (large 
grains) and the fluctuation dominated regime ( small grains) are qualitative different, but show that 
the crossover between them is very smooth with no signs of critical level spacings separating them. 
We compare our DMRG results with the exact ones obtained with the Richardson solution finding 
complete agreement. We also propose a simplified version of the DMRG wave function, called the 
Particle-Hole BCS ansatz, which agrees qualitatively with the DMRG solution and illustrates what 
is lacking in the PBCS wave function in order to describe correctly the crossover. Finally we present 
a new recursive method to compute norms and expectation values with the PBCS wave function. 

PACS number: 74.20.Fg, 74. 25. Ha, 74.80.Fp 



I) INTRODUCTION 

A fundamental question posed in 1959 by Anderson 
is "at what size of particles and what degree of scatter- 
ing will superconductivity actually cease" |jj . He argued 
that when the average level spacing d is of the order of 
the BCS gap A superconductivity must disappear. This 
old question was considered in the past by several au- 
thors (g|j3) and has been recently revived due to the ex- 
periments with ultrasmall Al grains performed by Ralph, 
Black and Thinkham (RBT) [Q. The experiments show 
the existence of a spectroscopic gap which can be driven 
to zero by application of magnetic fields. RBT also found 
a parity effect meaning that the magnitude of the spec- 
troscopic gap is larger for grains with an even number of 
electrons than for odd ones. 

From a theoretical point of view Anderson's question 
is challenging since it concerns the applicability of the 
standard BCS theory at nanometer scales jq]. Despite 
of some theoretical works using the grand canonical BCS 
wave function &RQ| it was soon realized that the descrip- 
tion of ultrasmall metallic grains calls for a canonical for- 
malism since the fluctuations in the electron number are 
strongly suppressed by charging effects [O-H . A canon- 
ical treatment of the BCS wave function has been known 
in Nuclear Physics for decades []15|-fl7|] ( for a review see 
IJbSl). The nucleus have a fixed number of fermions and 
the parity effects are clearly observable and interpreted 
theoretically. The ground state of the nucleus can be de- 
scribed by a wave function which is the projection of the 
BCS ansatz to a fixed number of fermions. This is the so 
called projected BCS ansatz (PBCS). The techniques for 
dealing with the PBCS wave function have been trans- 
lated to the study of ultrasmall metallic grains O] . The 



trouble with the BCS state and to a certain extent with 
the PBCS ansatz is that they are mean field approxima- 
tions which do not take care of the fluctuation effects that 
are supposed to be important for very small grains. An 
alternative is to use unbiased numerical methods where 
no assumption is made on the nature of the ground state. 
The authors of reference [|12| have studied systems up to 
25 electrons with the Lanczos method showing the im- 
portance of the logarithmic corrections in the supercon- 
ducting gaps first proposed in reference [lj] using a per- 
turbative renormalization group method. However exact 
diagonalization techniques cannot handled large systems 
where the crossover between the few-electron and the 
bulk superconducting regime is taking place for the ac- 
tual value of the BCS coupling constant, which for the Al 
grains is given approximately by A ~ 0.224 mU. Another 
alternative is to use the Density Matrix Renormaliza- 
tion Group Method (DMRG) Q which allows to study 
large systems with very high accuracy. This approach 
was initiated by the authors in reference PJJ, obtaining 
results which agree with those of the Lanczos method 
for small systems while improving the PBCS results for 
larger grains. In this paper we shall present a systematic 
study of the crossover region for grains with sizes in the 
range 20 up to 400, showing the importance of the fluc- 
tuations, which cannot be handle appropriately by the 
BCS or PBCS approaches. 

The BCS pairing Hamiltonian that we shall study in 
this paper has been solved exactly long time ago by 
Richardson in a series of papers between 1963 and 1977 
in the framework of Nuclear Physics |22j-|24[] . These pa- 
pers scaped the attention of the physics community until 
the recent developments in the field of ultrasmall metallic 
grains. Thus we have the great opportunity to compare 



the numerical results obtained with the DMRG method 
and the exact results obtained with the Richardson's 
wave function. Upon this comparison we shall see that 
the DMRG provides exact numerical results within a cer- 
tain accuracy which can be improved systematically by 
increasing the number of states kept. 

The overall picture we get from our study is that the 
few-electron and the bulk-limit regimes are qualitative 
different but the crossover is completely smooth. In 
this sense our results clarify and overcome the short- 
comings of previous grand-canonical BCS and canonical 
PBCS studies. In the BCS analysis superconductivity 
ceases to exists for level spacings d greater than a crit- 
ical value which is different for even grains cP c = 3.56A 
and for odd grains d\ = d° c /4 [|. In the PBCS study 
of Braun and von Delft the latter breakdown of super- 
conductivity does not occur but is replaced by a sharp 
crossover between the bulk regime and the fluctuation 
dominated regime which depends on the parity of the 
grains (d° c ~ 0.5A,d l c ~ 0.25 A) @. The results pre- 
sented in this paper will show no sign of critical level 
spacings separating qualitative different regimes. In fact 
we have been able to parametrize in a simple manner the 
numerical results found for several observables. These 
fitting-formulas are a sort of finite-size scaling similar to 
those that appear in low dimensional systems. 

The main tool we employ in our study is the particle- 
hole DMRG (PHDMRG) method first proposed in ref- 
erence 21 1 . This method follows the general philosophy 
of the real space DMRG method |2Q] but exploits the 
existence of a Fermi surface and the fluctuations around 
it. To apply the PHDMRG we have first to perform a 
particle-hole transformation where the Fermi sea is the 
vacuum of the basic operators. The states that appear in 
the DMRG are the particle-hole (p-h) excitations around 
the Fermi sea labelled by an integer £ that counts the 
number of particle pairs or holes pairs. Since we work 
at half filling, i.e. number of electrons equal to the num- 
ber of doubly degenerate states, the number £ is common 
to both particle and hole excitations in the ground state 
of the system. The DMRG algorithm selects the most 
probable p-h states that contribute to the exact ground 
state of the system. For every value of £ there are usually 
more than one p-h state, which form a sort of multiplet 
with multiplicity rag. The sum of all these multiplicities 
equals the total number m of states kept in the DMRG, 
i.e. m — Y]f rag. In our computations we have used a 
value of m = 60 which is sufficient to study system sizes 
up to 400 energy levels with a relative error of 10~ 4 in 
condensation energies. An outcome of the DMRG results 
is that for every value of £ there is a single p-h state which 
carries most of the probability. This fact suggests a sim- 
plified version of the DMRG based on an ansatz with only 
one p-h state per £. We call this state the particle-hole 
BCS ansatz (PHBCS). The reason for this terminology is 
that the PBCS state itself is a PHBCS state, though of a 



special type. While the PHBCS ansatz is a generic linear 
superposition of p-h states labelled by £, the PBCS state 
is a particular linear superposition of p-h states. We have 
thus a hierarchy of canonical variational ansatzs 



PBCS C PHBCS C DMRG C Exact 



(1) 



where every one contains its predecessor and is expected 
to give better results. From the PBCS to the PHBCS 
ansatzs one gains the freedom to mix different p-h states 
while in the DMRG ansatz, in addition to the latter free- 
dom, there are multiple p-h states for each value of £. 
We shall make a comparative analysis of the numerical 
results which will clearly show the qualitative and quanti- 
tative importance of these ingredients. The last member 
in the chain (P stands for the exact Richardson's solu- 
tion of the BCS model. We shall see that the numerical 
results obtained with the DMRG and the Richardson's 
solution are for practical purposes indistinguishable. 

The organization of the paper is as follows. In section 
II we define the model that is used to study ultrasmall 
metallic grains and summarize its essential features. In 
section III we introduce the PBCS wave function. In 
section IV we perform the p-h transformation, which is 
used to express the PBCS state in the p-h basis. We then 
propose the PHBCS state and find the effective Hamilto- 
nian that governs its dynamics. In section V we discuss 
in detail the DMRG method and relate it to the PHBCS 
ansatz. In section VI we present our numerical results for 
various quantities of interest obtained with the DMRG, 
PHBCS and PBCS methods. In section VII we state our 
conclusions. Technical details and derivations have been 
collected in two appendices. In appendix A we propose 
a novel recursion method to compute norms and expec- 
tation values with the PBCS state. In appendix B wc 
derive the form of the pairing BCS Hamiltonian in the 
p-h basis. 



II) THE BCS PAIRING HAMILTONIAN 

The BCS pairing Hamiltonian used for small metallic 
grains is given by |i| |lj] 

o n 

H = Yl fe ~ ^) c L c j> ~ Xd J2 c i+4,- c 3,- c j,+ 

j=l,<y=± »,J=1 

(2) 

where i,j = 1, 2, . . . , Q label single particle energy levels 
whose energies are given for simplicity by 6j — jd, where 
d is the average level spacing which is inversely propor- 
tional to the size of the grain. cj i<7 are electron destruc- 
tion operators of time reserved states a — ±. Finally 
fi is the chemical potential and A is the BCS coupling 
constant, whose appropriate value for the Al grains is 



0.224 JDJ. Given N e electrons they can form no Cooper 
pairs and b unpaired states such that N e = 2n + b. The 
number of electrons N e is equal to be number of states fi 
appearing in (0) . The Hamiltonian (0) decouples the un- 
paired electrons and hence & is a conserved quantity. The 
b unpaired electrons only contribute to the total ground 
state energy £ j, with their kinetic energy. Of particu- 
lar interest is the study of the parity effect which means 
that grains with an even number of electrons are more su- 
perconducting than odd grains. This phenomena, which 
occurs also in finite nuclei, can be characterized by the 
dependence of different observables as functions of b [|l8| . 
The Hamiltonian (|2|) has two regimes depending on 
the ratio d/A = 2sinh(l/A)/f2, between the level spacing 
d and the bulk superconducting gap A []6|— |l4|| . In the 
weak coupling region (d/A >> 1), which corresponds to 
small grains or small coupling constant, the system is in a 
regime with strong pairing fluctuations above the Fermi 
sea which lead to logarithmic renormalizations |y|. In 
the strong coupling regime (d/A << 1), which corre- 
sponds to large grains or strong coupling constant, the 
bulk-BCS wave function describes correctly the GS prop- 
erties. Using the grand canonical BCS wave function the 
crossover between the weak and strong coupling regimes 
occurs at d a c / A ~ 3.56 (even grains) and d\j 'A ~ 0.89 
(odd grains) [||. 



\PBCS(b)) = --i— re: o 6 +1 4,+ (itX Ivac) (7) 
ZnoMo = (vac|r' 2 X (ll ) ^ |vac) (9) 



While the PBCS state (0) depends on tf, variational 
parameters gi, the PBCS state (0) depends only on 2no 
parameters associated to the non-blocked levels. The un- 
paired states only contribute to the energy of the state 
(0) with the kinetic energy q. 

We can give a pictorial representation of the PBCS 
states (H) and (0) , which will be used later on in the dis- 
cussion of the DMRG. A system with non blocked levels, 
i.e. b = 0, can be represented as 



n/2 fi/2-i 



2 1 P 1 2 

• • o o 



fi/2-i a/2 
o o (10) 



where • denotes the p th particle level, o denotes the h th 
hole level and /x is the chemical potential separating par- 
ticles and holes. A system with one blocked level at the 
Fermi level is represented as 



n no — 1 



ft 



n -l n 
O O 



(11) 



III) THE PROJECTED BCS WAVE FUNCTION 

Let us first consider the case where all the electrons 
form Cooper pairs which can occupy all the allowed states 
of the system, i.e. $7 = 2no and 6 = 0. 

The PBCS wave function is given by 



\PBCS(b = 0)) 



n/2 



y/Zn, 



(4) Ivac) 



j n/2,n 



Vac|r n / 2 (rt) fV2 



vac) 



(3) 

(4) 
(5) 



where |vac) is the Fock vacuum of the electron operators 
and the variational parameters of the ansatz gi are re- 
lated to the standard BCS parameters Ui and Vi by the 
equation 



9i = — , 



v? = 1 



(6) 



The state (0) is the projection of the grand canoni- 
cal BCS state exp(r)|vac) into the Hubert space of Sl/2 
Cooper pairs. 

Let us consider now the case of b unpaired electrons. 
As explained in the previous section these electrons de- 
coupled from the rest of the system occupying the closest 
states to the Fermi level, namely i = no + 1, . . . ,no + b. 
The latter levels are also called blocked states. The 
PBCS state for b > is given by 



where no is the total number of Cooper pairs and ff is the 
unpaired spin lying on the Fermi level. Finally a system 
with 6 = 2 unpaired electrons will be represented as 



n n -l 



1r I li- 



no— 1 ™o 

O O 



(12) 



In what follows we shall concentrate on the case 6 = 0, 
leaving for the appendices the cases with 6 > 0. The 
variational parameters gi in the ansatz (|3j) and ([7|) are 
found by minimization of the mean value of the Hamil- 
tonian (|2]) . This requires the computation of the norm of 
the PBCS states and the expectation value of (||). This 
problem was first considered in Nuclear Physics where 
the projection of the BCS wave function was needed in 
order to take into account the finite size effects of the 
nucleus |[l5|Jl7|Jl^1 . The method developed in references 
|17| leads to a set of 2n coupled equations which are 
solved in terms of a set of auxiliary quantities entering 
the computation. In appendix A we propose an alterna- 
tive method based on recursion relations which can be 
easily implemented for system sizes f2 < 400. We have 
checked that this method reproduces the same results 
obtained by Braun and von Dclf |p~4| who used the tech- 
niques of references [ [17[ . The recursion method is quite 
manageable and will be used later on to study the PH- 
BCS ansatz. 



IV) THE PARTICLE-HOLE BCS STATE 

In the weak coupling limit d/A >> 1 the separation 
between energy levels is much greater than the bulk su- 
perconducting gap. The physics of this regime is given 
by the fluctuations around the Fermi state, 



\FS) 



a/2 



(13) 



= cj+ c l- ( see Appendix A for notations). 
An appropriate choice of the chemical potential /j, in (0) 
guarantees that particle and hole excitations around the 
Fermi sea (13) have the same energy. This symmetry 
implies that the PBCS parameters gt satisfy the following 
relation, 



gn+i- 



— , i = l,. 
9i 



,fi 



(14) 



which holds in particular for the BCS solution for the 
variational parameters Ui and Vi in eq. (ra). Eq.(|lJ) 
is a consequence of the particle-hole symmetry of the 
Hamiltonian (0) that we shall show more explicitly below. 



The PBCS state in the particle-hole basis 

In order to take full advantage of the symmetry condi- 
tion ( |l4]) it is convenient to establish the relationship be- 
tween the PBCS state (§) and the Fermi sea \FS). With 
this aim we shall write the pairing operator Tq given in 
©as 

T Ci =T A (x)+T B (l) (15) 

Fa(x) = 2J P =i x pPp 

r s(|) = Eft=l ^ P h 

where p,h = 1, . . . , $7/2 label the particle and holes states 
starting from the levels closest to the Fermi sea, i.e. 



Pp = Pn/2+p, Ph 



P 



fJ/2+l-hi 



(p,h= 1, 



,$7/2) 
(16) 



and x p = Xh(p = h) are the gi parameters for the particle 



states. 



9n/2-\ 



P=l, 



,$7/2 



(17) 



In ( p"5| ) we have used the eq.(|l4|). Eq. ( pif ) gives the 
transformation from the original pairing operators Pi to 
the new operators P p and Ph- While the vacuum state 
| vac) is annihilated by Pi, Vi, the Fermi state \FS) is 
annihilated by P p and P^ . Eq. ( |l6| ) is nothing but the p-h 
transformation used in BCS to go from the Fock vacuum 
to the Fermi sea. 



The operator T' A creates a pair of particles above the 
Fermi sea while the operator Tb creates a pair of holes. 
Hence we can use these operators to expand a basis of 
particle-holes states above the Fermi sea. Let us define 
the normalized state 



10 



1 



Zg,n/i{x) 



{T\{x)f (T B (x)Y\FS) 



(18) 



which is simply the tensor product of the particle state 
\£)a with I particles and the hole state \£)b with £ holes. 
One can show that the PBCS state (||) can be expanded 
in the p-h basis (|l8|) as follows p8|, 



\pbcs) = y.%1^ bcs \z) 



where 



t 



PBCS _ 



m/m 2 



v 7 ! 



S7/2 7 S7^f2/2,£7/2 



(x) 



ZlSl/2{x) 

(*!) a 



(19) 



(20) 



As a simple application of the formula (EOT) let us 
consider the PBCS state characterized by the choice 
x p = l,Vp, which corresponds to a fully superconducting 
state. The p-h amplitudes are given by 



t 



PBCS 



(Xi 



1) 



C, 



n/2,t 



N' 



a 



€l/2,Cl/2 



(21) 



where C NM = Nl/(M\(N - A/)!). This is an interest- 
ing result for it implies that the probability wi = ipj for 
finding the p-h state \£) in ^ PBCS is given by the hyper- 
geometric series distribution 



we 



r 2 

^Q./2,l 



C\ 



a,n/2 



a/2 

1=0 



1, (Xp = 1) 



(22) 



In the limit when $7 is large the distribution Q22 
becomes a normal distribution centered at $7/4 with 
quadratic deviation yT2/2. This result is the basis of 
the DMRG method applied in |||] to the pairing BCS 
Hamiltonian. 

Incidentally, it is interesting to observe that the distri- 
bution (E2|) is the same as the one found by Kaulke and 
Peschel for the S z = ground state of the Heisenberg 
ferromagnet pq j. The reason for this correspondence is 
based on the pseudo spin representation of the pairing 
Hamiltonian (g) ( see Appendix A). 

Eq.([19|) means that the PBCS state can be seen as the 
superposition of p-h states \£) with amplitudes tpf BCS , 
which both depend on the variational parameters x p . As 
explained in the introduction we can try to relax eq. (|19j) 
and consider ipi as variational parameters independent on 
the parameters x p . This will lead us to a more general 
ansatz which shares many common properties with the 
DMRG state. 



The Particle-Hole BCS Ansatz 

The previous study leads us to consider a general p-h 
state of the form. 



n/2, 



\phbcs) = Y:\Lli>t\e)A®\ti 



(23) 



The lowest energy of H G gives the ground state conden- 
sation energy divided by d. In appendix B we derive the 
Hamiltonian in the p-h basis for a general value of b. 

The p-h symmetry of the Hamiltonian (G9) amounts 
to its invariance under the following mappings, 



K- 



K B , A^ B 



(30) 



where \£) A and \£)b are the particle and hole pieces of the 
state given in eq. ( |l8| ) and ipt are independent parame- 
ters not constrained to satisfy eq.(|20|). Strictly speaking 
the p-h states ( p3| ) belong to the Hilbert space Hphbcs 
expanded by the p-h basis Pq) and their dynamics is 
governed by the projection of the pairing Hamiltonian 

In order to find this effective Hamiltonian acting in 
"Hphbcs it is convenient to express (0) using the p-h 
operators (pit), together with the p-h number operators, 



In the p-h basis \£) the Hamiltonian (|27j) becomes a 
tridiagonal matrix. This fact can be proved using the 
factorization of every state (23) into its particle and hole 
contents, 
given by 



The unique non-vanishing entries of H are 



{i 



(e\H c \e) 

-1\H C \£) 



2 A (£\(K A -XA^A)\£) A 
-X A (£-1\A\£) A 



(31) 



N„ 



2P p tp p , N h = 2P h Pl 



(24) 



A simple computation yields, 



The state if) a bas the same form as the PBCS state de- 
fined in eq.(g|) with the replacements gi — » x p , fl —> Q/2. 
Hence we can compute the matrix elements appearing 
in ([5l|) by using the auxiliary quantities introduced in 
Appendix A, 



H 



2J2 n h L\id( 



2 ' 
(f 



2 1 



l-ti)-n 

\-p)-l4N p 
+ EL=i [-d(% + l-h)+u + Xd] N h - 

HE py P$P P ' + Y,h,w PhPl + E Pjl (p^P h + PpPt)) 



(£\H C \£) 



(25) 



-MEp, pl (* P s e p, 



(l-l)xffi y 



(£-l\H c \£) = -X 



Zp-i 



[12pSp 



(32) 



This Hamiltonian has a p-h symmetry provided we 
choose the following chemical potential, 



/i = - (O + 1 - A) 



(26) 



which guarantees that the particle and hole excitations 
have the same energy. Using (p6|) the Hamiltonian (p£ 
adopts the simple form, 



H/d=-(§) 2 + K' 



K 1 



(27) 



A (AU + B^B + AB + A^B^) 



where 



K A = J2%1 e P N p , K B = ^Ll e h N h 
£p = £h=p-\ + \, (P = h) 
A = EpllP P , B = Y%L\Pl 



(28) 



The term — d (^) in ( p7| ) gives the energy of the Fermi 
sea with the chemical potential (f2q). We can subtract 
that term and measure the energy in units of d, 



H c = ('2 



H/d 



(29) 



The numerical procedure to find the PHBCS state with 
lowest energy is summarized in the following steps: 

i) Make an initial guess for the parameters x p . One 
can use for example the BCS values. 

ii) Construct the effective Hamiltonian (^2[) for this 
choice of parameters using the recursion method given in 
Appendix A. 

iii) Find the lowest GS of the effective Hamiltonian 

iv) Change slightly the parameters x p and repeat the 
steps ii) and iii), comparing the GS energy so obtained 
with the one determined in the previous step. Stop the 
process until convergence is achieved. 

Another important point is that in the PHBCS state 
defined in ( p3[ ) we can actually restrict the sum over £ 
to only a small number of values. For example we can 
include the states from up to say £ max and check the 
convergence in the energy by changing £ max . In the range 
il < 400 it is enough to choose £ max = 11. 

This method gives the values of x p and ipg of the PH- 
BCS state that minimizes the energy of the BCS pairing 
Hamiltonian. We shall present our results in section VI. 



V) THE DMRG STATE 

The DMRG state represents the next step in our 
route to go beyond the PBCS ansatz. Let us denote 



by {\a, &)a}™Li & n orthonormal set of me many body 
particle states containing £ particles, i.e. 



n+l 



y(a,£\a',£')A = St,# $ a ,a> 



(33) 



Similarly we shall introduce a set {\j3, £)bYrL\ °f many 
body hole states with £ holes. With these notations a 
DMRG state can be written as Cltl 



E E 1>aA*)\<*>*)A®\0>QB 



I Q,/3=l 



(34) 



Comparing eqs.(|23|) and ( p4| ) we see that the PHBCS 
states are a particular case of DMRG states where there 
is only one representative particle or hole state per £, 
namely 



i>. 



PHBCS 



<i™(4 K = 1,W) 



(35) 



A generic DMRG state involves higher multiplicities, 
i.e. me > 1, which is important for the numerical accu- 
racy of the method. 

Similar approximations to the DMRG in the context of 
strongly correlated systems have been given in references 

IH§- 

We shall next present the basic ideas of the DMRG 
method and its application to the pairing BCS Hamil- 
tonian pi]] . In the DMRG one has to break the system 
under study into two pieces called the system block A 
and the environment block B. In our case the block A 
contains all the particle levels while B contains the hole 
ones. If the system size, i.e. O, is large enough one can- 
not keep all the particle or hole states and hence one 
has to look for an effective description of them. This is 
done by keeping a set of m particle ( rcsp. hole) states 
\a,£)A, a= 1, ... ,me, \(3,£) B , P — 1, ■ ■ ■ ,me, as in eq. 
(pi) with, 



E 



m ( 



(36) 



These two sets of states are chosen in such a way that 
the state constructed in eq.(|34|) gives the best possible 
approximation to the exact GS of the whole system. The 
construction proceeds in successive steps starting from 
small grains. We begin with a system with = 4 en- 
ergy levels, which are chosen as the closest two parti- 
cle and hole states near the Fermi level /i. This system 
can be represented as • • oo, where we use the notation 
introduced in eqs. ( |10|Jll| , p^ ) . For larger systems, i.e. 
il = 2(n + 1) with n > 1, the whole system is described 
by the superblock •^4„i3 n o, where the block A n ( resp. 
B n ) gives an effective description of the n particle ( resp. 
hole) levels closer to the Fermi energy in terms of the m 
dimensional basis introduced above. In the notation of 
cqs.(PJll|,P) we have 





An 




n 


2 


1 


• 


• . • 


• 





s„ 




1 


2 


n 


o 


o ■ • 


• o 



n+l 
O 



(37) 



A generic state of the superblock •A n B n o 1 in the sector 
with equal number of particles and holes, reads 



I 1 / 1 ) = Ea,/3^ 8 Vw(4,4,4,4) X 

\h) n+x O \a,£ 2 ) An ® |/3,4)s n ® |4)n+i, 
(4+4 = 4 + 4) 



(38) 



where |4)n+i is the (n+l) th particle state which is empty 
for ^i = and occupied for £\ = 1. The hole state 
|4)n+i is similarly defined. The dynamics of the wave 
function (Bq) is governed by the superblock Hamiltonian 
which we shall construct below. The dimension of the 
Hilbert space of the superblock, diuiHsB, is smaller than 
4m 2 , for the constraint 4+4 = 4 + 4 eliminates many 
states. dmiHgB is usually much smaller than the exact 
dimension of the Hilbert space of states with £1 levels 
at half filling which is given by the combinatorial num- 
ber Cn,n/2- For example for il — 24 the latter number 
is 2704156, while the largest superblock matrix involved 
in the DMRG calculation with m = 60 has dimension 
3066. Another example is given by £1 = 400 where the 
dimension of the Hilbert space is of order 10 119 , while 
the largest superblock dimension is also 3066. 

The next step in the DMRG is to find the lowest 
eigenstate of the superblock Hamiltonian using the Lanc- 
zos technique. The corresponding eigenvalue gives the 
DMRG estimate of the GS energy for the system with 
£1 = 2(n + 1) energy levels. Since the DMRG is a varia- 
tional method it gives an upper bound of the exact result. 
Moreover the GS of the superblock previously found can 
be used to construct the new blocks A n +i and B n +i that 
give the effective description of the lowest n+l particle 
and hole states. This is achieved by first constructing the 
reduced density matrix of the subsystem »A n by tracing 
over the hole subsystem S n o, 

P'£, (44,fif 2 )= (39) 

E/VsA ^(4,4, 4, 4) ifia'Atu?*' 4,4) 



The density matrix (39) has a block diagonal form 
where each block is labeled by the total number of par- 
ticles, i.e. £ = £\ + £2- Let us denote the corresponding 
density matrix p* A . It is easy to see that it is a square 
matrix with dimension me + mi-\. One can also de- 
fine a reduced density matrix for the hole subsystem B n ° 
by tracing over the particle subsystem, however the p-h 
symmetry implies the equality of the particle and hole 
density matrices. This is a sort of reflection symmetry 
that recalls the symmetry between left and right blocks 
used in the infinite system DMRG algorithm applied to 
Id systems |2(|. In fact the particle-hole DMRG pro- 
posed above is an improved infinite system algorithm, 
obtained with some modifications to be explained below. 



Of course we can also deal with cases where the 
particle-hole symmetry does not hold. In this cases the 
particle and holes states kept in the DMRG will differ. 

Given the density matrix p' , we diagonalize it and 
find its eigenvalues 



( Wl {l) 



p\ A = O e 



w 2 {t) 



\ 



\ 



Of (40) 



J m£-\-m£ 



-A*) J 



where O is an orthogonal matrix and w\{l) > W2(£) > 
.... Once we have found all the eigenvalues for all al- 
lowed values of £ we put them together and sort them in 
decreasing order of magnitude. The DMRG truncation 
*A„ — > A' n +i consist in choosing the first m eigenvectors 
with highest eigenvalue. The renormalized block A' n +i 
will be described by a set of m' f states such 
that m = Y^,£ fn'i ( recall eq.(|36[)). The change of basis 
from the old block •An to the new block A 1 n +i is given by 
the first m'g column vectors of the orthogonal matrix Oe . 
The error of the truncation is measured by 1 — P m (P m = 

Em \ 

k=l w k)- 

Let us now give the Hamiltonian H.abo of the su- 
perblock •A n B n o, 



h»ab 

- Hab + H. A 



= H A + H B +H. + H t 



H Ac 



+ H.b + H Bo + H mo (41) 



H A 

H. 
Hab 

H.a 

H A o 

H.o 



K A - 

n 
^n+l^n+1 

—X(A n B n 



\A n A n 



h.c.) 



(p) f p (p) 



P. 



= -\{A n P^ + h.c.) 
= -\{A n P^ + h.c.) 

= -\(Pi%P n % f + h.c.) 



(42) 



where N, 



(p) 



Pi p) and P n h) are defined in eqs.@, © 



and (|24j). The superindices have been introduced to dis- 
tinguish between the particle and hole operators. The op- 
erators A n ,B n , K A and K n coincide with those defined 
in J2g| ) with 0/2 replaced by n. The terms Hb,H ,H,b 
and Hbo can be derived from those of (142) by the p- 
h transformation (pO]). The splitting ( |4l| ) of the su- 
perblock Hamiltonian H.abo recalls the one used by Xi- 
ang in the momentum space DMRG pit] and more re- 
cently by White and Martin in their study of the water 
molecule [B2j. However there are important differences 
between the latter approaches and ours. First of all Xi- 
ang's method uses a finite system algorithm while ours 
is an infinite system one combined with a renormaliza- 
tion of the interaction to be explained below. Secondly 
we exploit the p-h symmetry of the problem which is not 
the case of references |3l]|32| . 

The DMRG provides a many body description of the 
blocks A n and B n , which means that the operators acting 
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FIG. 1. GS condensation energy for fi = 100 and A = 0.4 
computed with the DMRG method as a function of the num- 
ber of states kept (i.e. m). The exact result is given by 
£(exact) = -40.5007557623. 

on these blocks are represented by m x m matrices. In 
our case the operators that we need to keep track are 
[A n ], [A n A n ] and [Nj]. The DMRG proposed above is 
an infinite system algorithm, which is sufficient to study 
moderate system sizes (N < 400). A way to improve the 
numerical accuracy of the infinite system method is to 
choose an effective value of the coupling constant A„ at 
the n th DMRG step in such a way that the value of the 
bulk gap is the one of the final system. This is guaranteed 
by the equation 



. . 1 2(n + 

smh A: ^ -n 



— sinh— 
A 



(43) 



VI) NUMERICAL RESULTS 

Comparison of the DMRG with exact results 

A system with $7 = 24 levels can be exactly diago- 
nalized with the Lanczos techniques as done in ref. [ [La . 
The DMRG calculation with m — 60 agrees with the 
exact Lanczos condensation energy in the first 7 digits. 
For larger systems the Lanczos method cannot be applied 
but as we said in the Introduction one can use the exact 
Richardson's solution. In fig.l we plot the exact GS con- 
densation energy for a system with £1 = 100 levels and 
A = 0.4, together with the DMRG results as a function 
of the number of states kept (i.e. m). One can clearly see 
the exponential convergence in m of the DMRG towards 
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FIG. 2. Condensation energies of the 6 = state as a func- 
tion of ft obtained with the DMRG, PHBCS and PBCS meth- 
ods. The energies are normalized respect to the bulk super- 
conducting gap given by A = d£7/(2sinh(l/A)). 

the exact solution. Another comparison we have made 
is for a system with H, = 400 and A = 0.224. Keep- 
ing m = 60 states we get for the GS condensation 
energy Eg (DMRG) /d = -22.5168 with an estimated 
relative error of 10~ 4 . The exact result is given by 
Eg (exact) /d = -22.5183141, which is within the esti- 
mated error. For lower system sizes the relative error is 
smaller so for all practical purposes the DMRG results 
cannot distinguished from the exact ones. In the figures 2 
through 8 presented below the curves labelled by DMRG 
also provide the exact results. 
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Q. 

FIG. 3. Condensation energies of the 6=1 state obtained 

with the DMRG, PHBCS and PBCS methods. 



• The DMRG method gives much lower condensation 
energies than those of the PBCS method, while the 
PHBCS lies in between ( see figs. 2 and 3). 

• The sharp crossover of the PBCS results, which is 
reflected in a sudden change in the slope of E c h for 
6 = and 1 as a function of O, is completely absent 
in the DMRG and the PHBCS results. 

• The dependence of Eg on il is rather smooth and 
can be parametrized by fitting the DMRG curves 
with the following formula ( see fig. 4 ) 



Eg /a = -a b n -p b + 7biog(n)/n 



(44) 



Condensation Energy 

The crossover between the superconducting and fluctu- 
ation dominated regimes can be neatly characterized by 
the condensation energy Eg defined as the difference be- 
tween the total energy £& of the GS and the energy of the 
uncorrelated Fermi sea l-FS 1 }- This energy has been com- 
puted for even and odd grains using the grand canonical 
(g.c.) BCS wave function JqJTn] and the canonical PBCS 
wave function |L4|. The g.c. studies suggest a break- 
down of superconductivity for large values of d while in 
the canonical case this breakdown is replaced by a sharp 
crossover between two different regimes at a characteris- 
tic level spacing dg ~ 0.5A. For d < dg the condensation 
energy Eg is an extensive quantity (~ 1/rf) correspond- 
ing to a BCS-like behaviour, while for d > dg the energy 
Eg is an intensive quantity (almost independent of d) 



In figures 2 and 3 we plot the DMRG, PBCS and PH- 
BCS results for the condensation energies Eg for even 
grains (6 = 0) with sizes ranging from 22 up to 400 
and odd grains (6=1) for sizes between 21 and 401. In 
figures 4 we collect the DMRG results corresponding to 
6 = 0,1,2 and 3. 

In these figures we observe the following features: 



where the constants a>b,(3b and 7b are given in ta- 
ble 1. The fitting formula fl44j) is an improved ver- 
sion to the one used in reference (21j and can be 
motivated from physical considerations as will be 
discussed below. 
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2 
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0.005701 
0.004586 
0.003439 

0.002747 


2.6678 
2.2463 
2.1258 
2.0485 


3.9321 
5.3275 
6.9290 
8.1536 



Table 1. Values of the parameters of formula ( [44[ ) 

that gives the best square least fit of the DMRG 

data plotted in fig. 4. 

The fit (H|) is specially good for the 6 = 
DMRG data but it is also quite performant for 
the other states 6 > 0. The first term in J44| ) 
represents the bulk correlation energy given by 
Eg = -A 2 /(2d),V6. Using the relation d/A = 
2sinh(l/A)/0 we deduce that the parameter a& 
should be independent of 6 taking the following 
value, 



1 



4sinh(l/A) 



= 0.005757 for A = 0.224 (45) 




FIG. 4. Condensation energies of the b = 0, 1, 2 and 3 states 
obtained with the DMRG method. The continuum lines are 
give by the fit ( |44[ ) with the numerical coefficients given in 
table I 

We see from table 1 that ao is close to the bulk 
value (|4q), while a;,, for b > 0, have not still reached 
that value. The constant term (3^ depends smoothly 

on b. This fact agrees with the computation of 
E^ using second order perturbation theory which 
yields 



[3 = 21n(2) A 2 sinh(l/A) = 3.0206 



(46) 



This value is close to those shown in table 1. The 
coefficient jb which controls the logarithmic term 
in (|4|) behaves roughly as 7^ — c\ + C26, where 
C\ = 3.9 and C2 = 1.4. This type of behaviour 
agrees qualitatively with second order perturbation 
theory, though the values of c\ and ci are different. 

In summary, eq.(Eih combines the extensive be- 
haviour (afe-term), the intensive behaviour (/3&- 
term) and the logarithmic corrections (7b-term) in 
a simple manner, showing no sign of sharp crossover 
in the condensation energy as a function of the 
grain's size. This conclusion is supported by fur- 
ther evidences shown below. 



Spectroscopic gaps: Parity effect 

The parity-dependent spectral gaps are defined as 



Eq = £ 2 — £0, (even grains) 
Ei = £3 — £\ , (odd grains) 



(47) 



In figure 5 we plot the DMRG, PHBCS and PBCS 
results, which we next comment. 

• All the results share the same qualitative features 
namely, E^ > Eq for 




FIG. 5. Spectroscopic gaps E® measured in units of d. 
The subscripts e and o correspond to the cases b = and 
6=1 respectively. 

Q, < fi c , while Ef < E§ for U, > fl c . The value of 
J7 C depends slightly on the method used, i.e. il c ~ 
200. 

• Quantitatively however the DMRG gives much 
greater spectroscopic gaps than the PBCS method, 
specially for the odd grains. 

• The difference Eq — Ef for Q, > 51 c is smaller for 
the DMRG than the PBCS method, which means 
that the parity effect is smoother in the former 
method. 



Matveev-Larkin's parameter 

Another characterization of the parity effect is in terms 
of a gap parameter which measures the difference be- 
tween the GS energy of an odd grain and the mean en- 
ergy of the neighbour even grains obtained by adding and 
removing one electron fll£fl , 

A ML - Slip.) - ~ (£b(fi + 1) + £o(^ - 1)) (48) 

In fig. 6 we display our results. 
Comments: 

• All the curves show a minimum in Aml/A as a 
function of d/A. This latter feature was first con- 
jectured by Matveev and Larkin JTy] and confirmed 
by Mastellone at al. JT3] using the Lanczos method 
and the PBCS method by Braun and von Delft pi • 

• The shape of the DMRG curve is rather smooth as 
compared with the PBCS and the PHBCS meth- 
ods. This can be interpreted as a suppression of the 
even-odd parity effect in agreement with the results 
found for the spectral gaps. 
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FIG. 6. Matveev-Larkin's parameter obtained with the 
DMRG, PHBCS and PBCS methods 



• The DMRG results of fig. 6 can be fitted with the 
following formula, which can be derived from the 
fits (p3) of the condensation energies, 



Aa/l/A = 0.4215 + 0.18375 f 
+ 0.09683 | - 0.01606 ^log^ 



(49) 



This eq. shows that in the region 0.3 < d/A < 3.5 
the logarithmic term is not very important. The 
logarithmic corrections are contained in the renor- 
malization of the coefficient of the term d/A, whose 
bare value is A/2 = 0.112. The constant term 
equals the difference (3 Q — (3\ of the condensation 
energies (see m4) and table 1). 



Pairing parameter 

The BCS superconducting order parameter is strictly 
zero in the canonical ensemble. For that reason one has 
to find another quantity to characterize the pair mixing 
across the Fermi level that takes place in the ground state 
for fixed number of electrons. We shall choose the pairing 
parameter proposed in references ]q,[l4]], 



A 6 = Ad£ i C,- 



(50) 



which measures the fluctuation in the occupation num- 
bers. In the g.c. BCS case Cj = UjVj and A& coincides 
with the usual superconducting parameter A. 

In figs. 7 and 8 we show our results for Ao and Ai 
respectively. 

Comments: 

• Fig. 7 shows that the sharp transition occurring in 
the PBCS ansatz between the strong and weak cou- 
pling regimes is completely absent in the DMRG 
state. In the latter state the pairing parameter, 
when measured in units of A, converges monoton- 
ically to its bulk limit from above. 
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FIG. 7. DMRG, PHBCS and PBCSresults for the Pairing 
parameter Aq as defined in equation 



• In the odd case the crossover predicted by the 
PBCS method is more dramatic than in the even 
one Q. The DMRG and PHBCS results show that 
this is an artifact of the PBCS ansatz. The exis- 
tence of a minimum for Ai and not for Ao is due 
to the blocking effect produced by the unoccupied 
single state at the Fermi level. 

• The PHBCS curves in figs 7 and 8 agree qualita- 
tively with the DMRG curves, while they differ 
strongly from the PBCS curves. This shows the 
importance of letting the amplitudes tpg to be in- 
dependent from the BCS-like parameters <7j. 



Particle-Hole Probabilities 

Another comparison between the DMRG, PHBCS and 
PBCS states can be given in terms of the probability of 
finding a state with £ particles or holes. If 4> is the GS 
of the whole system one has to construct the reduced 
density matrix for the particle or the hole subsystems 
and look for the corresponding eigenvalues. As shown 
in section IV the reduced particle density matrix of the 
PBCS and PHBCS states contains a unique eigenstate 
with probability W£ per number of particles £, given by 
wi = ipj, where ipt is given by eq. (EQ) for the PBCS 
state while ipt for the PHBCS has to be obtained through 
the minimization process explained at the end of section 
IV. In fig. 9 we display our numerical results for we as 
a function of fi. The reduced particle density matrix 
derived from the DMRG state has several eigenvectors for 
a fixed number of particles £, with eigenvalues w n (£)(n = 
1, . . .) ( see eq (|40|)). In fig. 10 we plot our numerical 
results for w n {£). 

Comments: 
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FIG. 8. DMRG, PHBCS and PBCS results for the Pairing 
parameter Ai as denned in equation ( pOJ ) 

• The overall pattern of the particle probabilities is 
common to all the ansatzs namely, i) the Fermi sea 
is the most probable state for < ft < il±, where 
the value of fli depends on the ansatz, ii) in the 
interval f2i < il < H.2 the most probable state has 
one particle, while the probability of the Fermi sea 
continue to decrease crossing over eventually the 
probability of a 2 particle state, iii) every curve 
associated to a given number of particles £, first 
increases for small grains, then reaches a maximum, 
where it is the most probable state, and then starts 
to decrease. 

• The probabilities of the PBCS states show the char- 
acteristic sharp crossover in the region 160 < 51 < 
220, in agreement with similar behaviour observed 
in the condensation energy Eq (fig.l) , spectro- 
scopic gap Eq (fig. 4) and pairing parameter A 
(fig.6). 

• In contrast to the latter behaviour, the PHBCS 
and DMRG probabilities evolve smoothly with the 
system size showing no signs of discontinuities or 
abruptness, in clear agreement with the observables 
computed above. 

• The PHBCS curves are in one to one correspon- 
dence with the most probable DMRG states, while 
the next most probable DMRG states, with the 
same number of particles, have much less proba- 
bility. This justifies a posteriori the PHBCS ansatz 
where multiple states with the same number of par- 
ticles are not included. 

• For a fixed system size the DMRG and the PHBCS 
probabilities decay roughly as wi ~ exp(— c\l — £ n \). 
This type of exponential decay has been observed 
also in DMRG studies of spins chain and explains 
the accuracy of the DMRG method since in that 
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a 

FIG. 9. Plot of the particle-hole probabilities we, = if)f for 
I = 0, 1, . . . , 7. The PHBCS ( resp.PBCS) results are given 
by the continnum (resp. discontinuum) curves. 

case a small number of states kept per block is 
enough to faithfully reconstruct the exact ground 
state. 

• Finally we observe in fig. 10 that the next most 
probable DMRG states reproduce essentially the 
same pattern as the most probable ones. The same 
is true for the next to next most probable ones and 
so on. There seems to be a sort of self-similar struc- 
ture whose origin would be interesting to under- 
stand. For fl very large we expect that all these 
states will have a very small probability so that 
only the most probable ones would be necessary to 
describe the GS. In this case the PBCS and PH- 
BCS should coincide asymptotically. To show that 
this happens we have to consider system sizes larger 
than those studied in this paper. 



VII) CONCLUSIONS 

The main conclusion we draw from the results pre- 
sented in the previous section is that the crossover be- 
tween the fluctuation dominated regime and the bulk 
limit is completely smooth in the sense that there are 
no critical level spacings separating a superconducting 
phase and a fluctuation dominated phase. This result 
clarifies and overcomes the short-comings of previous 
grand canonical and canonical BCS studies. The abrupt 
crossover obtained with the PBCS state is an artifact of 
that method. Our DMRG results agree with the exact 
solution with an accuracy of at least 10 -4 for condensa- 
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PB95-01123 (J.D.) and PB97-1190 (G.S.). 
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FIG. 10. Plot ol the particle-hole DMRG probabilities 

w n (t), which are defined as the eigenvalues of the reduced 
density matrix with I particles or holes. The thick continuum 
lines correspond to n = 1 and t — 0, 1, . . . , 9, the discontin- 
uum lines correspond to n — 2 and £ = 1, 2, . . . , 7 and the 
thin continuum lines correspond to n = 3 and I — 1, 2, . . . , 5. 

Instead of a breaking or suppression of superconductiv- 
ity for ultrasmall grains we rather observe that supercon- 
ductivity and fluctuations cannot be genuinely separated 
and that they gradually mix with the system size. 

We have explained in more detail the particle-hole 
DMRG proposed in reference |2l]] which can be applied 
not only to the reduced BCS Hamiltonian with arbitrary 
energy levels but also to Hamiltonians where the pair- 
ing coupling may be level dependent, i.e. A — -> Xij. In 
this sense we can in principle study, using the particle- 
hole DMRG, the effect of level statistics gj3^] and more 
general pairing interactions where no exact solution is 
available. 

We have developed a new recursive method to deal 
with the PBCS wave function which is somewhat simpler 
than the methods currently used. 

We have proposed a new wave function, the particle- 
hole BCS state (PHBCS), which stands somehow in be- 
tween the PBCS and DMRG states and which can be 
studied using the recursive method mentioned above. 
The PHBCS also shows a smooth crossover between large 
and small grains correctly describing the interplay be- 
tween superconducting correlations and fluctuations. 
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APPENDIX A) PBCS STATES : RECURSION 
RELATION METHOD 

In this appendix we shall present a new method to 
compute norms and expectation values of observables in 
the PBCS state (0). Let us first define the following 
operators, 



p} = 4+4 



Pi = p 



Nt = c 



+ c *,+ 



•ct a 



(51) 



which satisfy the commutation relations, 



p. pT 



= Sij (l - #«) 



Ni,P, 



= 26 Z1 P (52) 



Eqs.(p2[) imply that the pairing creation Pj, the pair- 
ing destruction Pi and the electron number JVj operators 
satisfy an SU(2) algebra. This is the basis of the pseudo 
spin representation of the Hamiltonian (0) which can be 
written as 



ff=x>-A*)#i-Adx; p i tp j 



(53) 



3=1 



i,j=l 



For the non-blocked levels we can make the replace- 
ments Pj — > af ,Pi — > a^ ,Ni — > (of + 1) and transform 
( p3| ) into a XY Hamiltonian with non local interactions 
and a position dependent magnetic field. 

The collective pair operator (f|) and condensate (|J) can 
be written as, 



4 = E^ p / . i^) = r 



iN 



vac) 



(54) 



In order to find the norm and the energy of the PBCS 
state (pj) we shall introduce the following auxiliary quan- 
tities 

Z N = ^ r t^ 



& - (rffpt r j*-i 



z» = (rff-^p/rir 1 



(55) 



iN 



^JV-2 



fJV 



PP^n 



where all the expectation values are computed respect 
to the vacuum state. Using the commutation relations 
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( p2| ) we derive the action of annihilation operators on 
the condensate 

N t \N) = 2N 9i P} \N - 1) (56) 

P, \N) = N 9l \N-1)-N{N-1) gfp} \N - 2) (57) 

The recurrence relations for the quantities defined in 
@ are 



Z&i =9i(N- 1) Sf" -gf(N- 1) (N - 2) T i} 



N-l 



APPENDIX B) THE PAIRING BCS 

HAMILTONIAN IN THE PARTICLE-HOLE 

BASIS 

In section IV we gave the expression of the Hamilto- 
nian (0) in the p-h basis. We shall derive below the corre- 
sponding expressions for arbitrary values of the blocked 
levels b. 

Using the operators (|I7j) and (J23) we can write the 
Hamiltonian (H) as, 







(58) 


z» = z N - 1 - 


-{N-\) gi Sf~ X 


(59) 


f = N 9l Z N - x - 


-N(N-l)g*S?- 1 


(60) 



"o+b "o 

H= J2 (ei-Ui) + 2j2(e h 

i— no + 1 h— 1 






(69) 



T$ = NgjSf - 1 -N(N-1) fiz"- 1 (61) 

The matrices Z and T are symmetric and T has null 
diagonal matrix elements. These properties are not ex- 
plicitly manifested in the recurrence relations ( pq-pl| ) . In 
order to make these properties evident we insert (p8[) and 
( p0| ) into (|6l]) obtaining 

9l g N(N-l)[Z N -^ (62) 



no no 



h=i 



-ME p l p v + E ^ + E (^ + p ^) ] 



p,p' 



h.h' 



ij 



(N - 2) faS? ~ 2 + 9j Sf- 2 ) + (N-2)(N- 3) g i9j T»~*] 
We now define the hated quantities 

sr 



where the particle hole energy levels are e p = d(no + b + 
p), en = d(no + 1 — h), with p, h = 1, . . . , no- The equal- 
ity between the particle and hole energies is achieved by 
choosing the chemical potential « as 



S? = 



Z 



N 



jiN _ ij 

ij — %N 



fi = d[n 



(63) 



6 + 1-A 



(70) 



in terms of which the energy of the normalized state <p4 
reads, 



in which case the Hamiltonian ( p9|) adopts the simple 
form 



E = 2NJ2 i (e i - rfgiS? - XdN £ y 9i Sf 
Equations (pOh and ( pi] ) are transformed into 

Z N 



-S 



N 



Ngi-N(N-l)gtS. 



2nN-l 



ZN-1 t 

^f^ = g i9j N{N-l){l- 



(64) 

(65) 
(66) 



H/d = -n (n + 6) + f + ^ + K * 
- A (AU + Bt£ + AS + A+B+) 



(71) 



where 



e P = e fc =p+*=|t*, (p=ft) (72) 



7V-2 



(N - 2) (^Sf - 2 + gjS»- J ) +(N-2)(N- 3) g^ 

Taking into account that '^2 li giS? = Z N , multiplying 
(p3) by gi and summing over i we get a relation for the 
norm ratios 

4^ = at E gf - n (n - 1) E gf sr 1 (67) 



The constant term in (j71l) gives the energy of the Fermi 
sea with the chemical potential (|70|). The correlation 
energy E^ in units of d is given by the lowest eigenvalue 



of the Hamiltonian H^ 



H£ = +H/d+ UoK + 6)-^ 



Z' 

Eqs. (H | 
Z° = l 



together with the initial conditions 



* = £< 



cl _ 9i 

^ ~Zl 



(68) 



(73) 



The total energy £{,(0) of a grain with fi electrons and 
b blocked levels can be obtained by adding the chemical 
potential term to (|7l|), 



r Q /n 



can be used to find the values of Sj and T/. that deter- 
mine the energy (p3) of the PHBCS state. 



2 V 2 



f6(n) = e^ (n) + d[- - + 1 - a + - - + a 



b (b 



2 V 2 



(74) 
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From this eq. we can easily relate the spectroscopic 
gaps E^ and the condensation energies E^ , 



E° = £ b+2 (n) - S b (Q) 



(75) 



Similarly the Matveev-Larkin gap parameter dehned 



in (48) can be obtained as, 

a ml = £i(n) - 1 (£o(n + 1) + £b(fi - 1)) (76) 
+ E c x {n) - \ (E c (n + 1) + E c (n - 1)) 



M 

2 
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